home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_8 / a8_3.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.5 KB  |  107 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 8.3 (Local Minimum Search: Quadratic Interpolation).
  9. % Section    8.1,    Minimization of a Function, Page 416
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program finds the minimum of a function f(x).
  13.  
  14. % The method is quadratic interpolation.
  15.  
  16. % It is assumed that f(x) is unimodal over [a,b].
  17.  
  18. % The function  f(x)  is to be placed in the M-file  f.m
  19. % function z = f(x)
  20. % z = x.^2 - sin(x);
  21.  
  22. delete f.m
  23. diary f.m; disp('function z = f(x)');...
  24.            disp('z = x.^2 - sin(x);');...
  25. diary off;
  26.  
  27. f(0); % Remark. f.m and quadmin.m are used for Algorithm 8.3
  28.  
  29. pause % Press any key to continue.
  30.  
  31. clc;
  32.  
  33. % Place the endpoints interval [a,b] in  a  and  b.
  34.  
  35. a  = 0;
  36.  
  37. b  = 1;
  38.  
  39. pause    % Press any key to graph the function.
  40.  
  41. clc; clg;
  42. hs = (b-a)/150;
  43. Xs = a:hs:b;
  44. Ys = f(Xs);  % User defined function.
  45. axis([-0.05 1.05 -0.25 0.16]);...
  46. plot(Xs,Ys,'g');...
  47. hold on;...
  48. plot([-0.05 1.05],[0 0],'b',[0 0],[-0.25 0.16],'b');...
  49. xlabel('x');...
  50. ylabel('y');...
  51. title('To search for the minimum of y = f(x).');...
  52. grid;...
  53. axis;...
  54. hold off;...
  55. shg; pause    % Press any key to continue.
  56.  
  57. clc;
  58.  
  59. % Place the endpoints interval [a,b] in  a  and  b.
  60.  
  61. % Place the tolerance for the abscissas in  delta.
  62.  
  63. % Place the tolerance for the ordinates in  epsilon.
  64.  
  65. a  = 0;
  66.  
  67. b  = 1;
  68.  
  69. delta = 1e-14;
  70.  
  71. epsilon = 1E-15;
  72.  
  73. [p,yp,dp,dy,P] = quadmin('f',a,b,delta,epsilon);
  74.  
  75. pause    % Press any key to find the minimum of f(x).
  76.  
  77. clc; clg;
  78. hs = (b-a)/150;
  79. Xs = a:hs:b;
  80. Ys = f(Xs);
  81. YP = f(P);
  82. axis([a b -0.25 0.175]);...
  83. plot(Xs,Ys,'g',P,YP,'or');...
  84. hold on;...
  85. plot([a b],[0 0],'b',[0 0],[-0.25 0.175],'b');...
  86. xlabel('x');...
  87. ylabel('y');...
  88. title('Searching for the minimum of y = f(x).');...
  89. grid;...
  90. axis;...
  91. hold off;...
  92. shg; pause    % Press any key to continue.
  93.  
  94. Mx1 = 'The results for a quadratic minimization search.';
  95. Mx2 = 'The solution is [p  f(p)];';
  96. Mx3 = 'The  abscissa  p ';
  97. Mx4 = 'error bound for p';
  98. Mx5 = 'The minimum value f(p)';
  99. Mx6 = 'error bound for f(p)';
  100. Mx7 = 'The list of iterations is:';
  101. Mx8 = '     p(k)               f(p(k))';
  102. clc,echo off,diary output,...
  103. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(''),...
  104. disp(Mx3),disp(p),disp(Mx4),disp(['± ',num2str(abs(dp))]),disp(''),...
  105. disp(Mx5),disp(yp),disp(Mx4),disp(['± ',num2str(abs(dy))]),...
  106. disp(''),disp(Mx7),disp(Mx8),disp([P;YP]'),diary off,echo on
  107.